Load preprocessed dataset (preprocessing code in data_preprocessing.Rmd)
glue('Number of genes: ', nrow(datExpr), '\n',
'Number of samples: ', ncol(datExpr), ' (', sum(datMeta$Diagnosis_=='ASD'), ' ASD, ',
sum(datMeta$Diagnosis_!='ASD'), ' controls)')
## Number of genes: 30109
## Number of samples: 80 (45 ASD, 35 controls)
5838 genes don’t have an adjusted p-value because they have less mean normalized counts than the optimal threshold link, so they are going to be considered not to be significant
plot_data = data.frame('id'=rownames(datExpr), 'mean_expression' = rowMeans(datExpr),
'SFARI_score'=DE_info$`gene-score`!='None',
'DExpressed'=DE_info$padj<0.05)
ggplotly(plot_data %>% ggplot(aes(x=mean_expression, fill=DExpressed, color=DExpressed)) + geom_density(alpha=0.3) +
scale_x_log10() + ggtitle('gene Mean Expression distribution') + theme_minimal())
We lose almost all of the genes with SFARI score
plot_data_SFARI = plot_data %>% filter(SFARI_score)
ggplotly(plot_data_SFARI %>% ggplot(aes(x=mean_expression, fill=DExpressed, color=DExpressed)) + geom_density(alpha=0.3) +
ggtitle('gene Mean Expression distribution for SFARI Genes') + scale_y_sqrt() + theme_minimal())
table(plot_data$DExpressed[plot_data$SFARI_score], useNA='ifany')
##
## FALSE TRUE <NA>
## 715 151 19
print(paste0('Losing ', round(100*(1-mean(plot_data$DExpressed[plot_data$SFARI_score==TRUE], na.rm=T)),1),
'% of the genes with SFARI score'))
## [1] "Losing 82.6% of the genes with SFARI score"
datExpr = datExpr[plot_data$DExpressed & !is.na(plot_data$DExpressed),]
DE_info = DE_info[plot_data$DExpressed & !is.na(plot_data$DExpressed),]
datGenes = datGenes[plot_data$DExpressed & !is.na(plot_data$DExpressed),]
rm(plot_data, plot_data_SFARI)
To make calculations more efficient for the more time consuming methods, we can perform PCA and keep the first 16 principal components which explain over 99.5% of the total variance
pca = prcomp(datExpr)
datExpr_redDim = pca$x %>% data.frame %>% dplyr::select(PC1:PC16)
clusterings = list()
clusterings[['SFARI_score']] = DE_info$`gene-score`
names(clusterings[['SFARI_score']]) = rownames(DE_info)
clusterings[['SFARI_bool']] = DE_info$`gene-score`!='None'
names(clusterings[['SFARI_bool']]) = rownames(DE_info)
clusterings[['syndromic']] = DE_info$syndromic==1
names(clusterings[['syndromic']]) = rownames(DE_info)
set.seed(123)
wss = sapply(1:10, function(k) kmeans(datExpr, k, iter.max=200, nstart=25,
algorithm='MacQueen')$tot.withinss)
plot(wss, type='b', main='K-Means Clustering')
best_k = 4
abline(v = best_k, col='blue')
datExpr_k_means = kmeans(datExpr, best_k, iter.max=100, nstart=25)
clusterings[['km']] = datExpr_k_means$cluster
Chose k=9 as best number of clusters. SFARI genes seem to group in the last two clusters
h_clusts = datExpr %>% dist %>% hclust
plot(h_clusts, hang = -1, cex = 0.6, labels=FALSE)
abline(h=19, col='blue')
best_k = 10
SFARI and Neuronal related genes seem to concentrate mainly in the aqua, blue and pink clusters
clusterings[['hc']] = cutree(h_clusts, best_k)
create_viridis_dict = function(){
min_score = clusterings[['SFARI_score']] %>% as.numeric %>% min(na.rm=TRUE)
max_score = clusterings[['SFARI_score']] %>% as.numeric %>% max(na.rm=TRUE)
viridis_score_cols = viridis(max_score - min_score + 1)
names(viridis_score_cols) = seq(min_score, max_score)
return(viridis_score_cols)
}
viridis_score_cols = create_viridis_dict()
dend_meta = DE_info[match(labels(h_clusts), DE_info$ID),] %>%
mutate('SFARI_score' = viridis_score_cols[`gene-score`], # Purple: 2, Yellow: 6
'SFARI_bool' = ifelse(`gene-score` != 'None', '#21908CFF', 'white'), # Acqua
'Syndromic' = ifelse(syndromic == T, 'orange', 'white'),
'Neuronal' = ifelse(ID %in% GO_neuronal$ID, '#666666','white')) %>%
dplyr::select(SFARI_score, SFARI_bool, Syndromic, Neuronal)
h_clusts %>% as.dendrogram %>% set('labels', rep('', nrow(datMeta))) %>%
set('branches_k_color', k=best_k) %>% plot
colored_bars(colors=dend_meta)
Samples are grouped into two big clusters, and then many outliers.
*The rest of the output plots can be found in the Data/Gandal/consensusClustering/genes_DE/ folder
Following this paper’s guidelines:
Run PCA and keep enough components to explain 60% of the variance (keeping 99.5% of the variance)
Run ICA with that same number of nbComp as principal components kept to then filter them
Select components with kurtosis > 3
Assign genes to clusters with FDR<0.01 using the fdrtool package
ICA_output = datExpr_redDim %>% runICA(nbComp=ncol(datExpr_redDim), method='JADE')
signals_w_kurtosis = ICA_output$S %>% data.frame %>% dplyr::select(names(apply(ICA_output$S, 2, kurtosis)>3))
ICA_clusters = apply(signals_w_kurtosis, 2, function(x) fdrtool(x, plot=F, verbose=F)$qval<0.01) %>% data.frame
ICA_clusters = ICA_clusters[colSums(ICA_clusters)>1]
clusterings[['ICA_min']] = rep(NA, nrow(ICA_clusters))
ICA_clusters_names = colnames(ICA_clusters)
for(c in ICA_clusters_names) clusterings[['ICA_min']][ICA_clusters[,c]] = c
clusterings[['ICA_NA']] = is.na(clusterings[['ICA_min']])
Leaves 76% of the genes without a cluster
ICA_clusters %>% rowSums %>% table
## .
## 0 1 2 3 4 5 6 7 8 9
## 2297 336 168 91 45 36 14 8 4 1
# ICA_clusters %>% mutate(cl_sum=rowSums(.)) %>% as.matrix %>% melt %>% ggplot(aes(Var2, Var1)) +
# geom_tile(aes(fill=value)) + xlab('Clusters') + ylab('Samples') +
# theme_minimal() + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) + coord_flip()
Trying again these time with all of the principal components and 40 clusters
ICA_output = pca$x %>% data.frame %>% runICA(nbComp=40, method='JADE')
signals_w_kurtosis = ICA_output$S %>% data.frame %>% dplyr::select(names(apply(ICA_output$S, 2, kurtosis)>3))
ICA_clusters = apply(signals_w_kurtosis, 2, function(x) fdrtool(x, plot=F, verbose=F)$qval<0.01) %>% data.frame
ICA_clusters = ICA_clusters[colSums(ICA_clusters)>1]
Doesn’t make a big difference (69%), but it’s still better
ICA_clusters %>% rowSums %>% table
## .
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 2074 283 188 135 95 63 49 43 30 19 6 8 4 3
clusterings[['ICA_min']] = rep(NA, nrow(ICA_clusters))
ICA_clusters_names = colnames(ICA_clusters)
for(c in ICA_clusters_names) clusterings[['ICA_min']][ICA_clusters[,c]] = c
clusterings[['ICA_NA']] = is.na(clusterings[['ICA_min']])
# ICA_clusters %>% mutate(cl_sum=rowSums(.)) %>% as.matrix %>% melt %>% ggplot(aes(Var2, Var1)) +
# geom_tile(aes(fill=value)) + xlab('Clusters') + ylab('Samples') +
# theme_minimal() + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) + coord_flip()
Note: This method does not work with the reduced version of datExpr.
best_power = datExpr %>% t %>% pickSoftThreshold(powerVector = seq(1, 30, by=2))
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.00342 0.261 0.988 7.75e+02 7.70e+02 1330.000
## 2 3 0.58100 -2.220 0.963 9.82e+01 9.03e+01 331.000
## 3 5 0.83600 -2.560 0.990 1.92e+01 1.58e+01 110.000
## 4 7 0.90000 -2.520 0.994 4.96e+00 3.49e+00 43.500
## 5 9 0.88700 -2.380 0.951 1.59e+00 9.14e-01 19.700
## 6 11 0.38500 -3.210 0.349 6.03e-01 2.66e-01 9.780
## 7 13 0.95700 -1.830 0.985 2.62e-01 8.50e-02 5.250
## 8 15 0.94200 -1.780 0.964 1.26e-01 2.95e-02 3.360
## 9 17 0.93400 -1.760 0.974 6.67e-02 1.08e-02 2.450
## 10 19 0.94800 -1.700 0.985 3.77e-02 4.11e-03 1.830
## 11 21 0.39700 -2.280 0.354 2.25e-02 1.65e-03 1.390
## 12 23 0.40700 -2.190 0.367 1.41e-02 6.55e-04 1.070
## 13 25 0.41200 -2.110 0.371 9.23e-03 2.72e-04 0.829
## 14 27 0.41400 -2.040 0.361 6.23e-03 1.15e-04 0.651
## 15 29 0.41300 -1.980 0.363 4.33e-03 4.94e-05 0.515
network = datExpr %>% t %>% blockwiseModules(power=best_power$powerEstimate, numericLabels=T)
clusterings[['WGCNA']] = network$colors
names(clusterings[['WGCNA']]) = rownames(datExpr)
It leaves 1211 genes without a cluster (~40%)
clusterings[['WGCNA']] %>% table
## .
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 1211 449 203 158 128 112 106 100 80 76 71 65 40 36 36
## 15 16 17 18 19
## 31 26 25 25 22
The BIC decreases monotonically, but it seems to stabilise at bit at 14
n_clust = datExpr %>% Optimal_Clusters_GMM(max_clusters=40, criterion='BIC', plot_data=FALSE)
plot(n_clust, type='l', main='Bayesian Information Criterion to choose number of clusters')
abline(v=14, col='blue')
best_k = 14
gmm = datExpr %>% GMM(best_k)
clusterings[['GMM']] = gmm$Log_likelihood %>% apply(1, which.max)
Separating the two clouds of points into two clusters
intercept=-0.2
slope=0.02
manual_clusters = as.factor(as.numeric(slope*datExpr_redDim$PC1 + intercept > datExpr_redDim$PC2))
names(manual_clusters) = rownames(datExpr_redDim)
clusterings[['Manual']] = manual_clusters
datExpr_redDim %>% ggplot(aes(PC1, PC2, color=manual_clusters)) + geom_point(alpha=0.3) +
xlab(paste0('PC1 (',round(100*summary(pca)$importance[2,1],1),'%)')) +
ylab(paste0('PC2 (',round(100*summary(pca)$importance[2,2],1),'%)')) +
geom_abline(intercept=intercept, slope=slope, color='gray') +
theme_minimal() + ggtitle('PCA')
clusterings[['Manual']] %>% table
## .
## 0 1
## 1382 1618
rm(intercept, slope, pca)
Both the aqua and the salmon clusters seem to be componsed of three Gaussians in the Mean and SD plots.
manual_clusters_data = cbind(apply(datExpr_redDim, 1, mean), apply(datExpr_redDim, 1, sd),
manual_clusters) %>% data.frame
colnames(manual_clusters_data) = c('mean','sd','cluster')
manual_clusters_data = manual_clusters_data %>% mutate('cluster'=as.factor(cluster))
p1 = manual_clusters_data %>% ggplot(aes(x=mean, color=cluster, fill=cluster)) +
geom_density(alpha=0.4) + theme_minimal()
p2 = manual_clusters_data %>% ggplot(aes(x=sd, color=cluster, fill=cluster)) +
geom_density(alpha=0.4) + theme_minimal()
grid.arrange(p1, p2, ncol=2)
Separate genes into four and two Gaussians, respectively by their mean expression:
gg_colour_hue = function(n) {
hues = seq(15, 375, length = n + 1)
hcl(h = hues, l = 65, c = 100)[1:n]
}
c1_mean = manual_clusters_data %>% filter(cluster==1) %>% dplyr::select(mean)
rownames(c1_mean) = rownames(manual_clusters_data)[manual_clusters_data$cluster=='1']
gmm_c1_mean = c1_mean %>% GMM(3)
c2_mean = manual_clusters_data %>% filter(cluster==2) %>% dplyr::select(mean)
rownames(c2_mean) = rownames(manual_clusters_data)[manual_clusters_data$cluster=='2']
gmm_c2_mean = c2_mean %>% GMM(3)
clusterings[['Manual_mean']] = as.character(clusterings[['Manual']])
clusterings[['Manual_mean']][clusterings[['Manual']]==0] = gmm_c1_mean$Log_likelihood %>%
apply(1, function(x) glue('1_',which.max(x)))
clusterings[['Manual_mean']][clusterings[['Manual']]==1] = gmm_c2_mean$Log_likelihood %>%
apply(1, function(x) glue('2_',which.max(x)))
plot_gaussians = manual_clusters_data %>% ggplot(aes(x=mean)) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(6)[1],
args=list(mean=gmm_c1_mean$centroids[1], sd=gmm_c1_mean$covariance_matrices[1])) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(6)[2],
args=list(mean=gmm_c1_mean$centroids[2], sd=gmm_c1_mean$covariance_matrices[2])) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(6)[3],
args=list(mean=gmm_c1_mean$centroids[3], sd=gmm_c1_mean$covariance_matrices[3])) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(6)[4],
args=list(mean=gmm_c2_mean$centroids[1], sd=gmm_c2_mean$covariance_matrices[1])) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(6)[5],
args=list(mean=gmm_c2_mean$centroids[2], sd=gmm_c2_mean$covariance_matrices[2])) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(6)[6],
args=list(mean=gmm_c2_mean$centroids[3], sd=gmm_c2_mean$covariance_matrices[3])) +
theme_minimal()
plot_points = datExpr_redDim %>% ggplot(aes_string(x='PC1', y='PC2', color=as.factor(clusterings[['Manual_mean']]))) +
geom_point() + theme_minimal()
grid.arrange(plot_gaussians, plot_points, ncol=2)
clusterings[['Manual_mean']] %>% table
## .
## 1_1 1_2 1_3 2_1 2_2 2_3
## 491 321 570 337 595 686
Separate clusters into three Gaussians per diagnosis by their sd:
n_clusters = 3
c1_sd = manual_clusters_data %>% filter(cluster==1) %>% dplyr::select(sd)
rownames(c1_sd) = rownames(manual_clusters_data)[manual_clusters_data$cluster=='1']
gmm_c1_sd = c1_sd %>% GMM(n_clusters)
c2_sd = manual_clusters_data %>% filter(cluster==2) %>% dplyr::select(sd)
rownames(c2_sd) = rownames(manual_clusters_data)[manual_clusters_data$cluster=='2']
gmm_c2_sd = c2_sd %>% GMM(n_clusters)
clusterings[['Manual_sd']] = as.character(clusterings[['Manual']])
clusterings[['Manual_sd']][clusterings[['Manual']]==0] = gmm_c1_sd$Log_likelihood %>%
apply(1, function(x) glue('1_',which.max(x)))
clusterings[['Manual_sd']][clusterings[['Manual']]==1] = gmm_c2_sd$Log_likelihood %>%
apply(1, function(x) glue('2_',which.max(x)))
plot_gaussians = manual_clusters_data %>% ggplot(aes(x=sd)) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(2*n_clusters)[1],
args=list(mean=gmm_c1_sd$centroids[1], sd=gmm_c1_sd$covariance_matrices[1])) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(2*n_clusters)[2],
args=list(mean=gmm_c1_sd$centroids[2], sd=gmm_c1_sd$covariance_matrices[2])) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(2*n_clusters)[3],
args=list(mean=gmm_c1_sd$centroids[3], sd=gmm_c1_sd$covariance_matrices[3])) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(2*n_clusters)[4],
args=list(mean=gmm_c2_sd$centroids[1], sd=gmm_c2_sd$covariance_matrices[1])) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(2*n_clusters)[5],
args=list(mean=gmm_c2_sd$centroids[2], sd=gmm_c2_sd$covariance_matrices[2])) +
stat_function(fun=dnorm, n=100, colour=gg_colour_hue(2*n_clusters)[6],
args=list(mean=gmm_c2_sd$centroids[3], sd=gmm_c2_sd$covariance_matrices[3])) +
theme_minimal()
plot_points = datExpr_redDim %>% ggplot(aes_string(x='PC1', y='PC2', color=as.factor(clusterings[['Manual_sd']]))) +
geom_point() + theme_minimal()
grid.arrange(plot_gaussians, plot_points, ncol=2)
clusterings[['Manual_sd']] %>% table
## .
## 1_1 1_2 1_3 2_1 2_2 2_3
## 461 337 584 452 603 563
rm(c1_sd, c2_sd, gmm_c1_sd, gmm_c2_sd)
# Clean up the environment a bit
rm(wss, datExpr_k_means, h_clusts, cc_output, best_k, ICA_output, ICA_clusters_names, signals_w_kurtosis,
n_clust, gmm, network, best_power, c, manual_clusters, manual_clusters_data, c1_mean, c2_mean,
gmm_c1_mean, gmm_c2_mean, p1, p2, dend_meta, plot_gaussians, plot_points, n_clusters, viridis_score_cols,
gg_colour_hue, create_viridis_dict)
Using Adjusted Rand Index:
Clusterings are not very similar except for K-Means, Hierarchical clustering and Gaussian Mixtures Model (which all do vertical clusterings)
No clustering method resembles the SFARI scores at all
Manual + Mean and Manual + SD resemble a bit K-means clustering and Hierarchical clustering
cluster_sim = data.frame(matrix(nrow = length(clusterings), ncol = length(clusterings)))
for(i in 1:(length(clusterings))){
cluster1 = as.factor(clusterings[[i]])
for(j in (i):length(clusterings)){
cluster2 = as.factor(clusterings[[j]])
cluster_sim[i,j] = adj.rand.index(cluster1, cluster2)
}
}
colnames(cluster_sim) = names(clusterings)
rownames(cluster_sim) = colnames(cluster_sim)
cluster_sim = cluster_sim %>% as.matrix %>% round(2)
heatmap.2(x = cluster_sim, Rowv = FALSE, Colv = FALSE, dendrogram = 'none',
cellnote = cluster_sim, notecol = 'black', trace = 'none', key = FALSE,
cexRow = 1, cexCol = 1, margins = c(7,7))
rm(i, j, cluster1, cluster2, cluster_sim)
The simple clustering methods only consider the 1st component, dividing by vertical lines
GMM does vertical clusters when using the complete expression matrix but round, small clusters when using the reduced version
WGCNA clusters don’t seem to have a strong relation with the first principal components
SFARI genes seem to be everywhere (perhaps a bit more concentrated on the right side of the plot)
1st PC seems to reflect the average level of expression of the genes
There seems to be a change in behaviour around PC1=0 (CC)
plot_points = datExpr_redDim %>% data.frame() %>% dplyr::select(PC1:PC3) %>%
mutate(ID = rownames(.), k_means = as.factor(clusterings[['km']]),
hc = as.factor(clusterings[['hc']]), cc = as.factor(clusterings[['cc_l1']]),
ica = as.factor(clusterings[['ICA_min']]),
n_ica = as.factor(rowSums(ICA_clusters)), gmm = as.factor(clusterings[['GMM']]),
wgcna = as.factor(clusterings[['WGCNA']]), manual = as.factor(clusterings[['Manual']]),
manual_mean = as.factor(clusterings[['Manual_mean']]), manual_sd = as.factor(clusterings[['Manual_sd']]),
SFARI = as.factor(clusterings[['SFARI_score']]), SFARI_bool = as.factor(clusterings[['SFARI_bool']]),
syndromic = as.factor(clusterings[['syndromic']])) %>%
bind_cols(DE_info[DE_info$ID %in% rownames(datExpr_redDim),]) %>%
mutate(avg_expr = log2(rowMeans(datExpr)+1)[rownames(datExpr) %in% rownames(datExpr_redDim)])
rownames(plot_points) = plot_points$ID
selectable_scatter_plot(plot_points, plot_points)
clusterings_file = './../Data/Gandal/clusterings.csv'
if(file.exists(clusterings_file)){
df = read.csv(clusterings_file, row.names=1)
if(!all(rownames(df) == rownames(datExpr))) stop('Gene ordering does not match the one in clusterings.csv!')
for(clustering in names(clusterings)){
df = df %>% mutate(!!clustering := sub(0, NA, clusterings[[clustering]]))
rownames(df) = rownames(datExpr)
}
} else {
df = clusterings %>% unlist %>% matrix(nrow=length(clusterings), byrow=T) %>% t %>% data.frame %>% na_if(0)
colnames(df) = names(clusterings)
rownames(df) = rownames(datExpr)
}
write.csv(df, file=clusterings_file)
rm(clusterings_file, df, clustering)
sessionInfo()
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## Running under: Scientific Linux 7.6 (Nitrogen)
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
##
## locale:
## [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
## [5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
## [7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid stats4 parallel stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] dendextend_1.10.0 gplots_3.0.1
## [3] pdfCluster_1.0-3 WGCNA_1.66
## [5] fastcluster_1.1.25 dynamicTreeCut_1.63-1
## [7] ClusterR_1.1.8 fdrtool_1.2.15
## [9] moments_0.14 MineICA_1.22.0
## [11] fastICA_1.2-1 Hmisc_4.1-1
## [13] Formula_1.2-3 survival_2.43-3
## [15] lattice_0.20-38 annotate_1.60.1
## [17] XML_3.98-1.11 Rgraphviz_2.26.0
## [19] igraph_1.2.4 colorspace_1.4-1
## [21] mclust_5.4 marray_1.60.0
## [23] limma_3.38.3 cluster_2.0.7-1
## [25] GOstats_2.48.0 graph_1.60.0
## [27] Category_2.48.1 Matrix_1.2-15
## [29] AnnotationDbi_1.44.0 IRanges_2.16.0
## [31] S4Vectors_0.20.1 gtools_3.5.0
## [33] biomaRt_2.38.0 xtable_1.8-3
## [35] foreach_1.4.4 scales_1.0.0
## [37] plyr_1.8.4 Biobase_2.42.0
## [39] BiocGenerics_0.28.0 JADE_2.0-1
## [41] ConsensusClusterPlus_1.46.0 GGally_1.4.0
## [43] gridExtra_2.3 viridis_0.5.1
## [45] viridisLite_0.3.0 RColorBrewer_1.1-2
## [47] plotlyutils_0.0.0.9000 plotly_4.8.0
## [49] glue_1.3.1 reshape2_1.4.3
## [51] forcats_0.3.0 stringr_1.4.0
## [53] dplyr_0.8.0.1 purrr_0.3.1
## [55] readr_1.3.1 tidyr_0.8.3
## [57] tibble_2.1.1 ggplot2_3.1.0
## [59] tidyverse_1.2.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.5 robust_0.4-18
## [3] RSQLite_2.1.1 htmlwidgets_1.3
## [5] trimcluster_0.1-2.1 BiocParallel_1.16.6
## [7] gmp_0.5-13.5 munsell_0.5.0
## [9] codetools_0.2-15 preprocessCore_1.44.0
## [11] withr_2.1.2 knitr_1.22
## [13] rstudioapi_0.7 geometry_0.4.0
## [15] robustbase_0.93-0 labeling_0.3
## [17] FD_1.0-12 GenomeInfoDbData_1.2.0
## [19] bit64_0.9-7 generics_0.0.2
## [21] xfun_0.5 diptest_0.75-7
## [23] R6_2.4.0 doParallel_1.0.14
## [25] GenomeInfoDb_1.18.2 clue_0.3-57
## [27] locfit_1.5-9.1 flexmix_2.3-15
## [29] bitops_1.0-6 reshape_0.8.7
## [31] DelayedArray_0.8.0 assertthat_0.2.1
## [33] promises_1.0.1 nnet_7.3-12
## [35] gtable_0.2.0 Cairo_1.5-9
## [37] rlang_0.3.2 genefilter_1.64.0
## [39] splines_3.5.2 lazyeval_0.2.2
## [41] acepack_1.4.1 impute_1.56.0
## [43] broom_0.5.1 checkmate_1.8.5
## [45] yaml_2.2.0 abind_1.4-5
## [47] modelr_0.1.4 crosstalk_1.0.0
## [49] backports_1.1.2 httpuv_1.5.0
## [51] RBGL_1.58.2 tools_3.5.2
## [53] Rcpp_1.0.1 base64enc_0.1-3
## [55] progress_1.2.0 zlibbioc_1.28.0
## [57] RCurl_1.95-4.10 prettyunits_1.0.2
## [59] rpart_4.1-13 SummarizedExperiment_1.12.0
## [61] haven_1.1.1 magrittr_1.5
## [63] data.table_1.12.0 mvtnorm_1.0-7
## [65] whisker_0.3-2 matrixStats_0.54.0
## [67] mime_0.6 hms_0.4.2
## [69] evaluate_0.13 readxl_1.1.0
## [71] compiler_3.5.2 KernSmooth_2.23-15
## [73] crayon_1.3.4 htmltools_0.3.6
## [75] later_0.8.0 mgcv_1.8-26
## [77] pcaPP_1.9-73 geneplotter_1.60.0
## [79] rrcov_1.4-3 lubridate_1.7.4
## [81] DBI_1.0.0 magic_1.5-9
## [83] MASS_7.3-51.1 fpc_2.1-11.1
## [85] ade4_1.7-11 permute_0.9-4
## [87] cli_1.1.0 gdata_2.18.0
## [89] GenomicRanges_1.34.0 pkgconfig_2.0.2
## [91] fit.models_0.5-14 foreign_0.8-71
## [93] xml2_1.2.0 XVector_0.22.0
## [95] AnnotationForge_1.24.0 rvest_0.3.2
## [97] digest_0.6.18 vegan_2.5-2
## [99] rmarkdown_1.12 cellranger_1.1.0
## [101] htmlTable_1.11.2 GSEABase_1.44.0
## [103] kernlab_0.9-27 shiny_1.2.0
## [105] modeltools_0.2-22 nlme_3.1-137
## [107] jsonlite_1.6 pillar_1.3.1
## [109] httr_1.4.0 DEoptimR_1.0-8
## [111] GO.db_3.7.0 prabclus_2.2-7
## [113] iterators_1.0.9 bit_1.1-14
## [115] class_7.3-14 stringi_1.4.3
## [117] blob_1.1.1 DESeq2_1.22.2
## [119] latticeExtra_0.6-28 caTools_1.17.1
## [121] memoise_1.1.0 ape_5.3